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ABSTRACT 

The  process  of  initial  orbit  determination,  or  catalogue  maintenance,  using  a  set  of  unlabelled  observations 
requires  a  method  of  choosing  which  observation  was  due  to  which  object.  Realities  of  imperfect  sensors 
mean  that  the  association  must  be  made  in  the  presence  of  missed  detections,  false  alarms  and  previously 
undetected  objects.  Data  association  is  not  only  essential  to  processing  observations,  it  can  also  be  one  of  the 
most  significant  computational  bottlenecks. 

The  constrained  admissible  region  multiple  hypothesis  filter  (CAR-MHF)  is  an  algorithm  for  initial  orbit 
determination  using  short-arc,  optical  (angles  only),  observations  of  space  objects.  CAR-MHF  uses  joint 
probabilistic  data  association  (JPDA),  a  well-established  approach  to  multi-target  data  association.  A  recent 
development  in  the  target  tracking  literature  is  the  use  of  graphical  models  to  formulate  data  association 
problems.  Using  an  approximate  inference  algorithm,  belief  propagation  (BP),  on  the  graphical  model  results 
in  an  algorithm  that  is  both  computationally  efficient  and  accurate. 

This  paper  compares  association  performance  on  a  set  of  deep-space  objects  with  CAR-MHF  using  JPDA 
and  BP.  The  results  of  the  analysis  show  that  by  using  the  BP  algorithm  there  are  significant  gains  in 
computational  load,  with  negligible  loss  in  accuracy  in  the  calculation  of  association  probabilities. 

1.  INTRODUCTION 

Data  association  is  an  essential,  but  challenging,  component  of  all  multiple  object  tracking  algorithms.  Where 
an  object  may  not  be  seen  (due  to  the  characteristics  of  the  object,  the  sensor  or  the  environment),  where 
false  measurements  are  produced  (due  to  a  noisy  sensor  or  other  environmental  effects)  or  where  there  are 
multiple  objects  or  multiple  measurements,  there  will  be  ambiguity  over  which  measurements  belong  to  which 
objects.  This  association  ambiguity  needs  to  be  resolved  in  order  to  incorporate  the  measurement  information 
into  estimates  of  the  object  states. 

The  constrained  admissible  region  multiple  hypothesis  filter  (CAR-MHF)  [1,  2,  3]  is  an  algorithm  for 
initial  orbit  determination  and  catalogue  maintenance  which  uses  the  joint  probabilistic  data  association 
(JPDA)  algorithm  [4],  JPDA  is  a  multi-target  data  association  algorithm  that  generates  an  association 
probability  for  each  measurement  with  each  object.  This  is  achieved  through  exhaustively  evaluating  the 
probability  of  all  valid  object-to-measurement  combinations.  The  version  of  JPDA  presented  in  this  paper 
additionally  accounts  for  the  mixture  of  Gaussians  describing  the  probability  density  of  the  state  of  each 
object  within  CAR-MHF. 

In  cases  where  there  are  many  objects  and  many  measurements  the  JPDA  algorithm  can  be  computationally 
expensive  [5].  There  are  practical  methods,  such  as  gating,  which  can  help  to  reduce  the  size  of  the  problems 
[4].  However  for  many  tracking  systems,  including  CAR-MHF,  JPDA  remains  a  computational  bottleneck. 
Approximation  methods  are  required  for  large  problems,  where  the  time  needed  to  explicitly  calculate  JPDA 
is  impractical.  A  promising  technique  for  approximating  the  measurement-to-object  association  probabilities 
is  the  use  of  belief  propagation  (BP)  within  a  graphical  model  of  the  association  problem  [6,  7,  8].  In  difficult 
air-defence  problems  BP  has  been  shown  to  be  much  faster  than  JPDA,  while  maintaining  sufficient  accuracy. 
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This  paper  aims  to  extend  the  analysis  of  the  BP  approximation  algorithm  to  orbit  determination  problems 
within  CAR-MHF. 

The  structure  of  this  paper  is  as  follows:  Section  2  discusses  the  constrained  admissible  region  (CAR) 
and  filtering  with  multiple  range  and  range-rate  hypotheses  (MHF).  Section  3  outlines  the  joint  probabilistic 
data  association  (JPDA)  algorithm  within  the  multiple  hypothesis  framework.  The  belief  propagation  (BP) 
approximation  to  JPDA  is  presented  in  Section  4.  Indicative  results  from  a  small  simulation  are  shown  in 
Section  5  with  conclusions  given  in  Section  6. 

2.  CAR-MHF 

2.1  The  Constrained  Admissible  Region 

The  constrained  admissible  region  (CAR)  is  an  area  in  measured  range  and  range-rate  found  by  constraining 
both  the  energy  (essentially  the  duration)  and  the  eccentricity  (the  shape)  of  possible  orbits  [1,  2].  The  CAR, 
as  discussed  in  this  paper,  relies  on  the  measurement  of  angles  and  their  rates.  Augmenting  the  angles  and 
angle-rates  with  feasible  ranges  and  range-rates  completes  the  six-dimensional  state  required  for  initial  orbit 
determination  giving  a  set  of  possible  positions  and  velocities  of  the  orbiting  object. 

The  energy,  £.  of  an  orbit  can  be  written  in  terms  of  the  inertial  position,  r,  and  velocity,  r,  of  the  object 
through 


where  /.t  is  the  standard  gravitational  parameter  of  Earth.  This  equation  can  be  rearranged  to  form  a 
relationship  between  the  measured  range  and  range-rate.  The  resulting  equation  [2]  produces  contours  of 
constant  orbital  energy  as  a  function  of  range  and  range-rate.  Similarly,  eccentricity,  e,  is  related  to  angular 
momentum,  h,  and  specific  energy  through 


Again  the  equation  can  be  rearranged  to  form  a  relationship  between  the  measured  range  and  range-rate 
resulting  in  contours  of  constant  eccentricity  as  a  function  of  range  and  range-rate. 

An  example  of  the  two  types  of  constraints  is  shown  in  Figure  1.  The  plot  on  the  left  of  Figure  1  is  the 
CAR  of  a  simulation  from  a  fictitious  sensor  located  at  DSTO  Edinburgh  observing  Optus  Cl  on  11th  July 
2014  at  15:45UTC.  The  plot  on  the  right  in  Figure  1  shows  a  more  representative  CAR  for  the  purposes  of 
initial  orbit  determination,  where  the  eccentricity  and  semi-major  axis  have  been  constrained  to  orbits  in 
near-GEO  (GEosynchronous  Orbit). 

2.2  Multiple-Hypothesis  Filter 

The  CAR  provides  a  constrained  region  of  the  object  state  space  which  is  supported  by  a  single  measurement. 
As  further  measurements  are  received  the  object  state  can  be  refined  with  that  information.  A  principled 
way  of  approaching  this  is  through  recursive  state-space  estimation.  For  this  application  a  non-linear  filter  is 
required  due  to  the  non-linear  nature  of  orbital  propagation  and  the  non-linearity  of  the  relationship  between 
the  measurements  and  the  orbital  state. 

The  unscented  Kalman  filter  (UKF)  [9]  is  used  in  this  application  to  propagate  the  probability  distribution 
of  the  object  over  time  and  refine  it  with  measurement  information.  While  the  UKF  is  typically  a  computa¬ 
tionally  efficient  and  accurate  filter  it  assumes  that  the  probability  of  the  state  is  Gaussian  distributed.  In 
general  a  CAR,  such  as  that  shown  on  the  right  in  Figure  1,  cannot  be  accurately  represented  by  a  Gaussian. 
To  ameliorate  this  shortcoming  a  mixture  of  multiple  Gaussians  is  used  to  represent  a  uniform  distribution 
over  the  CAR  to  any  desired  accuracy  [2].  This  mixture  results  in  a  filter  which  is  a  bank  of  many  UKFs,  one 
for  each  hypothesis  (component  of  the  mixture).  Clever  hypothesis  management  (which  is  not  reviewed  in 
this  paper)  results  in  the  number  of  elements  in  the  mixture  reducing  over  time  as  more  measurements  refine 
the  location  of  the  object. 

Following  the  derivation  of  the  well-known  Gaussian  sum  filters  [10,  11],  a  multiple-hypothesis  filter  can 
be  derived  which  accommodates  the  Gaussian  mixture  representation  of  the  CAR.  In  this  paper  Xk  is  used 


Figure  1:  The  constrained  admissible  region  (CAR).  The  figure  on  the  left  shows  the  CAR  for  several  values 
(Earth  radii)  of  the  semi-major  axis,  a,  in  blue  and  eccentricity,  e  in  red.  The  figure  on  the  right  shows  the 
CAR  with  semi-major  axis  between  42.15  x  106m  and  42.2  x  106m  (blue  lines)  and  eccentricity  between 
zero  and  ICR3  (red  ellipse).  The  grey-shaded  area  is  the  resulting  admissible  region  and  the  true  range  and 
range-rate  of  the  object  is  shown  by  the  cross. 


to  represent  the  state  of  the  object  at  time  fc,  zk  is  the  measurement  at  time  k  and  Zk  is  the  history  of  all 
measurements  Zk  =  {zi, . . . ,  zk}-  The  filter  follows  these  steps 

Initialisation  The  probability  distribution  of  the  object  at  time  A:  —  1  is  given  by  the  weighted  Gaussian 
mixture 

P(*fc_1|zfc-1)  =  E4-1^-1;4-nfc-1I^-1|fc-1),  (3) 

i 

where  A/"(a;;m,  E)  represents  a  Gaussian  distribution  in  variable  x  with  mean  to  and  covariance  E.  The 
factors  alk_1  represent  the  weight  on  each  component  of  the  mixture  defined  such  that  ak-i  =  1- 

Prediction  The  predicted  mixture  can  be  written  as 

pfalZ1*-1)  =Y^°‘k-iN{xk;xllk_1,Pllk_1)  (4) 

i 

where  xlk\k_1  and  PLk_1  can  be  calculated  using  the  standard  unscented  filter  formulation,  the  details 
of  which  can  be  found  in  [9].  Note  that  the  mixture  weights  are  unchanged  in  the  prediction  step. 


Update  The  mixture,  updated  with  a  measurement  at  time  k,  zk,  can  be  written  as 

p(xk\Zk)  =  Y,aW(xk;Zilk,Pilk) 


(5) 


where,  xk,k  and  Phfe  can  be  calculated  from  the  ith  predicted  component  and  the  measurement  at  time 
k  using  the  unscented  filtering  equations.  The  updated  component  weight  is  given  by 

dl.  =  alk_xA f{zk;zlk,Slk) 

i  ® k 


E 


j  ‘ 


where  zk  and  Sk  are  the  predicted  measurement  and  innovations  covariance. 


(6) 

(7) 


3.  JOINT  PROBABILISTIC  DATA  ASSOCIATION  IN  CAR-MHF 


The  MHF  presented  in  the  previous  section  assumes  that  there  is  a  single  object  and  a  perfect  sensor, 
that  is  no  false  alarms  and  no  missed  detections.  In  a  realistic  scenario  an  algorithm  must  be  able  to  deal 
with  multiple  objects,  false  alarms  and  missed  detections.  A  data  association  algorithm,  which  associates 
measurements  with  objects,  is  used  to  deal  with  these  practicalities.  In  this  case  we  use  the  joint  probabilistic 
data  association  algorithm  (JPDA)  [4]. 

The  imperfections  of  the  sensor  and  it’s  detection  performance  are  expressed  through  the  probability  of 
detection,  Pd  and  the  spatial  density  of  false  alarms,  A.  The  probability  of  detection  is  modelled  through 
the  characterisation  of  the  sensor  and  may  take  into  account  properties  of  the  object  (although  that  detail 
is  not  included  in  this  derivation  of  the  algorithm).  The  false  alarm  rate  is  typically  chosen  to  reflect  the 
performance  of  the  sensor,  but  is  also  related  to  the  processing  chain  that  produces  detections.  In  this  case 
the  number  of  false  alarms  is  assumed  to  be  Poisson-distributed. 


3.1  Multiple  object  MHF  with  association  ambiguity 

In  the  remainder  of  the  paper  Nk  is  the  number  of  objects  at  time  k  and  x\.  is  the  state  of  object  t  at  time 
k.  There  are  m measurements  at  time  k  and  so  now  the  history  of  all  measurements  is  Zk  =  {Z\, . . . ,  Zk} 
where  Zk  =  {zj., . . . ,  z™k}.  The  probability  distribution  of  the  state  of  object  t  is  represented  by 


p(xl\Zk)  =J2a)?Af(xtk;x)?lk,Pl 


i,t  \ 
k\k) 


(8) 


where  is  the  number  of  mixture  components  representing  the  probability  distribution  of  object  t. 

Define  [i'lf'  as  the  probability  that  measurement  j  g  {0,1,...,  m*,}  is  from  object  t,  where  j  =  0  corresponds 
to  the  object  not  being  detected.  Given  we  can  modify  the  update  step  of  the  MHF  in  the  previous 
section  to  reflect  the  ambiguous  measurement  origin.  The  mixture  updated  with  a  set  of  measurements  at 
time  k  can  be  written  as  (8)  where 


j.i.t  i.t 

°i  =  ak- 1  x 


PdX  zkf,  S'{’4)  for  j  >  0  (object  detected) 

(1  —  Pd)  for  j  =  0  (object  not  detected) 


a 


j.t  \  ^  j.i.t 

K  =  2^  a  k 

i= 1 

mk  oj,t 

;*  =  £  %<•'* 
3=  0  Wk 


(9) 

(10) 

(11) 


In  analogy  with  the  previous  section  z\f  and  Slk‘  are  the  predicted  measurement  and  innovations  covariance 
for  mixture  component  i  of  object  t.  Here  QjZ’4  is  the  *th  component  weight  for  the  jth  measurement  from 
the  tth  object,  vjf  is  the  likelihood  of  measurement  j,  given  that  it  came  from  object  t.  The  updated  weights 
are  a  sum  over  contributions  from  each  measurement  according  to  the  association  probabilities  /3^’4. 

In  a  similar  way,  and  Pk^k  can  be  calculated  as  a  weighted  sum  of  contributions  from  each  measurement. 
An  unscented  filter  gives  the  mean,  x’k)k  and  covariance,  PjjkX  of  the  updated  distribution  for  each 
measurement  from  each  component  of  the  mixture  from  each  object 


p(4’V 


Zk~')=N(x\ 


i,t  t 
k  ’ 


ck\k » 


rk\k  )• 


(12) 


In  the  special  j  =  0  (not  detected)  case  the  mean  and  covariance  come  from  the  predicted  distribution, 
^°k\k  ~  *fc|fc_i  an<^  =  Pk\k-v  Then  the  component  distribution  can  be  found  by  summing  over 

contributions  from  each  measurement 


p(a 


3=0 


i.t  ~i,i,t 
Jk  '■’Xk\k  ’ 


1  k\k  k 


(13) 


where 


(14) 


ft*  < 


W 


ah 


By  moment-matching,  we  can  calculate  the  updated  mean  and  covariance  of  each  component  in  an  object’s 
mixture  [4] 


i,t  _ 

k\k 


mk 


E 


w 


k  xk\k 


(15) 


3=0 

mk 

pk\k  =  2 L  < 

3=0 

Note  that  with  the  association  probabilities  /?£’*  supplied,  the  updated  distributions  above  are  calculated 
independently  for  each  object.  All  dependencies  that  arise  due  to  the  association  ambiguity  between  objects 
and  measurements  are  captured  in  the  association  probabilities.  As  such  the  computational  complexity  of  the 
algorithm  outlined  in  this  section  is  limited  to  0(ntkrrik)  for  each  object. 


p3,i,t 
r k\k 


[xk\k  ~xh\k 


)( 


j,i,t  ~i,t  ' 
Xk\k  ~  Xk\k, 


(16) 


3.2  Calculation  of  the  association  probabilities 

There  are  several  methods  that  could  be  chosen  to  calculate  the  association  probabilities  /3£f  [8].  CAR-MHF 
currently  uses  the  joint  probabilistic  data  association  (JPDA)  algorithm  [4],  which  is  outlined  here. 

For  each  object  t  €  {1, . . . ,  Nk}  we  define  an  association  variable  ak  €  {0, 1, . . . ,  m*,},  the  value  of  which 
is  an  index  to  the  measurement  with  which  the  object  is  hypothesised  to  be  associated  (zero  if  the  object  is 
hypothesised  to  have  not  been  detected).  Then  we  can  define  the  joint  association  variable  over  all  Nk  objects 

Ak  =  {al,...,a^ }.  (17) 

JPDA  proceeds  by  calculating  the  probabilities 

P{Ak\Zk)  (18) 


for  all  feasible  associations  Ak.  The  definition  of  a*,  ensures  that  for  a  given  Ak  each  object  can  be  associated 
with  a  single  measurement  (or  zero) .  A  feasible  association  must  also  be  constrained  such  that  a  measurement 
can  have  only  one  source,  that  is 

Nk 

Vj  >  0,  (19) 

i=l 


where  6(a,  b)  =  1  if  a  =  b  and  0  otherwise. 

Following  [4]  (with  some  rearrangement),  the  posterior  probability  of  a  joint  association  event  can  be 
written  as 

P{Ak\Zk}  (20) 

t=  1 


where  is  given  by 


lb3,t=  P^ftizAk)} 

k  \1-Pd 


for  j  >  0 
for  j  =  0, 


(21) 


and  ft  [^-(A;)]  is  the  likelihood  of  measurement  j  given  that  it  came  from  object  t.  Note  the  stark  similarity 
between  in  (21)  and  the  weights  akl,t  in  (9).  In  fact  for  the  case  of  a  Gaussian  mixture  ij)  should  be 
replaced  by  w  from  (10),  giving  the  probability  of  the  joint  association  event  [12] 


Mk 

P{Ak\Zk}  <x\[way . 

t= 1 


(22) 


Having  calculated  the  probability  of  all  of  the  joint  events,  P{Ak\Zk},  the  association  probabilities  can  be 
found  as  the  marginals  of  the  joint  distribution.  The  marginal  probability  of  an  association  event  is  given  by 


P{a[\Zk}=  PWl--.,a^\Zk}, 

and  then 

Pit  =  P{atk=j\Zk}. 


(23) 


(24) 


4.  BELIEF  PROPAGATION 


It  is  apparent  from  the  previous  section  that  the  computationally  complex  component  of  JPDA  is  the 
calculation  of  the  marginal  association  probabilities  through  the  explicit  evaluation  of  the  joint  association 
events.  JPDA  is  closely  related  to  the  #P-complete  problem  of  calculating  the  matrix  permanent  [5].  Gating 
[4],  which  involves  removing  unlikely  association  hypotheses  before  computing  association  probabilities,  is 
one  effective  measure  to  reduce  the  number  of  joint  events,  but  cannot  completely  remove  the  issue.  For 
large  problems  explicit  calculation  will  be  impossible  within  practical  time  and  memory  constraints,  so 
approximations  are  required.  This  section  presents  an  approximation  to  JPDA  based  on  belief  propagation 
(BP)  within  a  directed  graph.  The  BP  algorithm  was  developed  in  [6],  where  it  was  shown  that  this 
approximation  is  both  fast  and  accurate. 

The  previous  section  defined  the  association  variables  a\,  which  are  an  index  to  the  measurement 
hypothesised  for  object  t.  The  associations  can  be  equivalently  described  through  b°k  €  {0, 1, ... ,  Nk}.  for 
each  measurement  j,  the  value  of  which  is  an  index  to  the  object  with  which  the  measurement  is  hypothesised 
to  be  associated  (zero  if  the  measurement  is  hypothesised  to  be  either  a  false  alarm  or  a  new  object).  The 
sets  of  association  variables  are  redundant,  that  is  they  exhaustively  encode  the  association  hypotheses  in 
complementary  ways.  Rather  than  specify  the  constraints  on  the  association  variables  ak,  given  in  (19),  the 
BP  algorithm  considers  both  sets  of  association  variables  and  forces  consistency  between  them. 

The  JPDA  algorithm  calculates  the  marginal  probabilities,  P{atk\Zk},  explicitly  as  outlined  in  the 
previous  section.  The  belief  propagation  algorithm  approximates  these  probabilities  as  marginals  of  the  joint 
distribution 


P{al . . . ,  af ,  bl, . . . ,  b™*\Zk}  =  n  {  n  K) 


(25) 


The  factors 


al  =  3,  K  ±  t  or  b{  =  t,  a\  ±  j 
otherwise 


(26) 


ensure  the  consistency  between  the  a\  and  the  b?k  sets  of  association  variables.  The  ipkk,t  are  as  defined  in 
the  previous  section  and,  as  before,  should  be  replaced  with  the  relevant  sum  over  terms  in  the  Gaussian 
mixture.  The  desired  marginals  are  given  by 


P{a{\Zk} 


P{c 


'k’)  ' 


,iVfc 


bl,...,b™*\zk}. 


(27) 


The  graphical  model  used  to  encode  the  association  variables  is  shown  in  Figure  2.  The  belief  propagation 
algorithm  generates  messages  which  are  passed  between  nodes  along  the  edges  in  the  graph.  In  the  case 
of  a  tree-structured  graph  BP  is  exact,  however  for  a  graph  which  admits  cycles,  such  as  Figure  2,  BP  is 
necessarily  an  approximation.  The  message  passing  alternates  between  two  sets  of  messages;  those  passed 
from  the  left  to  the  right  (referring  to  Figure  2),  which  are  written  as  (u,  ,ty! )  then  those  passed  from  the 
right  to  the  left,  (%>  .  t ). 

uk  uk 


Figure  2:  The  graphical  model  used  to  encode  the  association  variables. 


As  shown  in  [7]  the  message  update  equations  simplify  to 

= 


1  Pd,  +  V’fc  Ab'-J-t 

1 


Ab->t  — 


i  +  E  oVi'^yj 


(28) 

(29) 


where  the  notation  (fit^yj)  has  been  used  for  (/u  t  . j  ).  Convergence  of  the  algorithm  is  guaranteed,  with  a 
proof  presented  in  [7],  and  a  numerical  test  for  convergence  of  the  algorithm  is  discussed  in  [8].  Once  the 
algorithm  has  converged  the  approximate  marginal  probabilities  are  given  by 


k  Ab-*-* 


fc  Ab'->* 


(30) 


5.  SIMULATION 


This  section  considers  a  cluster  of  26  objects  in  GEO  orbit.  The  CAR-MHF  framework,  using  the  JPDA 
algorithm  of  Section  3,  is  run  on  a  sequence  of  observations  of  these  objects  over  several  days.  In  the 
observation  instants  where  there  are  two  or  more  objects  competing  for  measurements,  the  performance  of 
JPDA  and  BP  is  compared.  A  total  of  192  association  examples  result  from  this  simulation.  The  comparison 
of  the  two  algorithms  is  made  through  both  computational  cost  and  accuracy  of  the  association  probabilities, 
/3^’4.  In  the  following  results  the  computational  cost  has  been  normalised  by  the  lowest  execution  time  of  the 
JPDA  algorithm. 

In  this  simulation  CAR-MHF  uses  a  false  alarm  intensity,  A  =  100,  and  the  probability  of  detection, 
Pd  =  0.95.  The  number  of  false  measurements  is  of  the  order  of  1  or  2  per  scan,  with  the  relatively  high  false 
alarm  intensity  resulting  from  the  precision  of  the  measurements.  It  should  be  noted  that  low  false  alarm 
rate  and  high  probability  of  detection  conditions  are  where  the  BP  algorithm  performance  is  weakest  [7].  As 
Pd  increases  and  A  decreases  the  BP  algorithm  decreases  in  accuracy  and  increases  in  compute  time. 

Figure  3  shows  the  histograms  of  the  (normalised)  computation  cost  for  both  JPDA  and  BP  algorithms 
on  the  same  192  association  examples.  It  is  immediately  apparent  that  the  BP  algorithm  is  not  only  quicker 
than  the  JPDA  algorithm,  but  its  run-time  is  more  consistent.  The  total  run-time  of  the  JPDA  is  2870,  while 
the  total  run-time  of  the  BP  algorithm  is  130,  indicating  that  BP  is  roughly  twenty  times  faster,  on  average. 
While  Figure  3  shows  the  distributions  over  compute  times,  it  does  not  indicate  the  relative  speed  of  the 
algorithms  on  specific  examples.  Figure  4  shows  the  ratio  of  the  JPDA  compute  time,  Tj,  to  the  BP  compute 
time,  Tb,  for  each  example,  r  =  Tj/Tb ,  so  values  r  >  1  imply  that  BP  is  faster,  while  r  <  1  imply  BP  is 
slower.  It  is  apparent  that  for  the  vast  majority  of  the  examples  BP  is  significantly  faster  than  JPDA,  but 
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Figure  3:  Histograms  comparing  the  compute  time  (normalised  to  the  shortest  JPDA  compute  time)  of  the 
JPDA  and  BP  algorithms. 
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Figure  4:  Histogram  of  the  ratio  of  JPDA  to  BP  compute  times.  Values  greater  than  1  imply  that  BP  is 
faster. 
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Figure  5:  The  magnitude  of  the  maximum  difference  between  the  BP  approximate  association  probabilities 
and  JPDA  in  the  non-negligible  cases  (6  out  of  192  examples). 


there  is  a  single  example  where  the  compute  time  for  BP  is  slower  (as  5  times)  than  explicit  JPDA  calculation. 
This  slow  example  also  corresponds  to  the  least  accurate  BP  approximation. 

Here  we  use  a  similar  metric  to  [7]  to  assess  the  accuracy  of  the  association  probabilities.  For  each  example 
we  find  the  maximum  difference  between  the  BP  calculated  marginal  probabilities  and  those  calculated  by 
JPDA.  In  the  majority  of  cases  the  differences  are  very  small,  less  than  0.01,  which  we  assume  to  be  negligible. 
There  are,  however,  six  cases  (about  3%  of  the  total)  which  have  larger  errors.  Each  of  these  errors  is  shown 
in  Figure  5.  The  largest  error  is  a  difference  in  probability  of  0.41.  This  corresponds  to  the  difference  between 
an  object  being  assigned  roughly  evenly  to  two  measurements  (JPDA)  and  the  object  being  assigned  to  a 
single  measurement  (BP).  It  is  the  subject  of  ongoing  analysis  to  determine  the  effect  of  these  differences  on 
overall  track  accuracy. 

In  all  but  one  of  the  cases  with  high  error  the  problem  is  small  (3  or  4  objects  assigned  to  2  or  3 
measurements)  and  the  computational  expense  in  calculating  these  examples  exactly  is  similar  to  BP.  However 
in  the  largest  case  (13  objects  being  assigned  to  7  measurements)  the  error  is  just  0.13  while  BP  is  roughly 
100  times  faster  than  JPDA.  Although  there  are  too  few  instances  in  this  simulation  to  make  a  conclusive 
judgement,  these  results  suggest  that  a  pragmatic  approach  would  be  to  apply  JPDA  to  small  association 
problems,  where  there  is  little  compuational  advantage  in  BP,  and  apply  a  BP  approximation  to  large 
association  problems  where  there  can  be  significant  computational  gains. 

Recent  research  [13]  has  shown  that  the  fractional  free  energy  proposed  in  [14]  can  provide  approximate 
marginal  probabilities  with  improved  accuracy  in  comparison  to  BP  for  problems  with  high  probability  of 
detection  and  low  false  alarm  rate.  Preliminary  work  has  shown  that  the  worst-case  error  can  be  reduced 
from  0.41  to  0.002  using  this  method  (below  our  threshold  for  negligible  difference),  although  the  method  in 
[13]  is  somewhat  slower  than  BP  for  the  problem  of  interest.  Rapid  solution  of  problems  with  this  formulation 
is  a  topic  of  ongoing  research. 

6.  CONCLUSIONS  AND  FURTHER  WORK 

This  paper  has  given  a  brief  introduction  to  the  constrained  admissible  region  multiple  hypothesis  filter 
(CAR-MHF)  for  initialising  and  maintaining  orbital  estimates  from  angle  and  angle-rate  measurements.  Joint 
probabilistic  data  association  (JPDA)  has  been  described  as  an  approach  to  deal  with  measurement-to-object 
association  ambiguity  within  CAR-MHF.  As  JPDA  can  be  a  computational  bottleneck  in  CAR-MHF,  a 
method  of  approximating  the  JPDA  calculations  has  been  evaluated,  based  on  belief  propagation  (BP)  in  a 
graphical  model  (see  Figure  2).  Comparisons  of  JPDA  and  BP  methods  on  a  simulated  data  set  have  shown 
that  in  the  majority  of  cases  BP  offers  a  significant  advantage  in  computation  time  (roughly  20  times  faster 
on  average)  with  negligible  difference  in  the  resulting  association  probabilities.  The  results  suggest  that  the 
small  number  of  cases  where  BP  produced  a  higher  error  could  be  largely  avoided  by  using  explicit  JPDA  for 
small  problems  and  limiting  the  use  of  BP  to  large  problems,  where  it  displays  the  most  benefit. 

We  have  two  ongoing  areas  of  research  related  to  data  association  within  CAR-MHF.  The  first  of  these  is  to 
analyse  the  effect  of  the  BP  approximations  on  the  track  quality  of  larger  simulated  data  sets.  While  admitting 
that  BP  is  an  approximation,  track  quality  is  the  ultimate  measure  of  system  performance,  regardless  of  the 
data  association  scheme.  The  second  area  is  in  the  investigation  of  alternative  schemes  for  approximating 
association  probabilities.  If  our  track  quality  analysis  concludes  that  the  BP  approximations  are  insufficient 
then  there  are  other  options  that  can  be  explored,  such  as  the  previously  mentioned  fractional  free  energy 
approach  [13]. 
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